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An  orthotropic  material  model  was  presented  and  its  applications  were  studied  in  this  report.  A  four-parameter  peridynamic 
model  for  orthotropic  materials  was  proposed  to  coincide  with  the  four  independent  material  properties  required  for  orthotropic 
materials.  This  four-parameter  model  was  different  from  other  models  as  it  had  four  parameters  and  was  mesh  independent. 
The  model  was  verified  by  a  static  tensile  test  and  a  vibration  excitation  of  a  laminated  beam.  An  SEN  (single  edge  notch)  test 
of  a  0°  laminated  plate  was  simulated  by  peridynamics  and  the  computational  results  matched  very  well  with  published 
experimental  results.  Fracture  initiation  and  crack  path  of  laminated  plates  with  different  fiber  orientations  were  also  studied 
using  peridynamics.  The  mesh-free  peridynamic  model  was  convenient  and  efficient  since  there  was  no  need  to  have  different 
meshes  for  different  fiber  orientations.  Without  using  special  meshes  and  requiring  prior  knowledge  of  fracture  paths,  a  general 
peridynamic  code  was  able  to  predict  fracture  velocity  and  crack  path  successfully. 

Technology  Transfer 


Peridynamic  Applications  for  Orthotropic  Materials 


It  has  been  shown  earlier  that  peridynamic  simulation  of  damage  process  does  not  require  any 
knowledge  of  the  damage  location  and  orientation  prior  to  the  simulation.  This  is  fundamentally 
different  from  finite  element  analysis  which  requires  knowledge  of  damage  location  and 
orientation  in  advance  to  impose  special  finite  element  mesh,  such  as  initial  damage  elements 
and  cohesive  zone  layers  [1],  for  damage  simulations.  This  prerequisite  becomes  even  more 
challenging  when  inhomogeneous  and  anisotropic  composite  materials  are  of  interest.  In 
addition,  peridynamic  simulation  does  not  require  remeshing  at  the  end  of  each  damage 
processing  step  since  it  is  a  mesh  free  method.  On  the  contrary,  finite  element  analysis  does. 
Based  on  these  difference,  peridynamics  should  be  more  suitable  for  simulating  dynamic  damage 
process  in  composite  materials  which  have  different  properties  in  different  locations  and  different 
orientations. 

Quite  some  simulations  of  composite  damage  process  have  been  available  in  the  literatures. 
Dwivedi  [1]  modeled  the  propagation  of  single-edge  notch  (SEN)  in  0°  laminated  plate  using 
cohesizve  zone  method.  Xu  [2]  and  Hu  [3]  proposed  a  two-parameter  discrete  peridynamic 
model  for  composite  damage  simulations,  in  which  there  were  two  kinds  of  bonds:  fiber  bond 
and  matrix  bond.  Two  material  properties,  a{  and  a2,  were  associated  with  the  two  types  of 
bonds.  Only  the  bonds  along  the  fiber  direction  were  associated  with  the  material  property  a, 
while  all  other  bonds  with  the  material  property  a2 .  This  model  required  remeshing  for  different 
fiber  directions.  For  example,  a  0°-90°  grid  mesh  could  only  be  used  for  a  0°  or  90°  laminae.  For 
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a  45°  lamina,  a  grid  mesh  consisting  of  45°  and  135°  was  required.  The  two-parameter  model 
was  an  aproximation  of  the  four  material  properties  invovled  in  ortho  tropic  materials.  They  were 
mainly  associated  with  two  Young’s  moduli,  E±  and  E2.  Its  capability  of  modeling  shear  behavior 
is  unknown. 

In  this  study,  a  continuous  orthrotropic  material  model  is  proposed.  It  is  based  on  continuous 
trigonometric  functions.  With  the  continuous  material  property  functions,  it  is  not  necessary  to 
have  bond  in  fiber  direction  and  therefore,  this  model  is  mesh  independent. 

5,1  Bar  model  for  orthotropic  materials 

This  model  is  based  on  the  bar  model  presented  in  Seciton  4.1.  The  peridynamic  equation  of 
motion  [4]  in  two-dimensional  domain  can  be  expressed  as 

pix  —  ffdA  +  b  (5.1) 

where  b  is  external  force.  The  force  boundary  condition  can  be  included  in  the  external  force. 

For  bar  model,  the  bond  function  /  is 

f  =  c-s  (5.2) 

where  s  is  bond  stretch.  Contrary  to  the  isotropic  material  model,  bond  material  property  c  is 
assumed  to  be  a  trigonometric  function 

c  —  d±  cos(0  —  a )4  +  d2  co s(0  —  a)2  +  d3  (5.3) 

where  d1}  d2  and  d3  are  constants  and  can  be  identified  from  composite  material  properties.  6  is 
bond  direction  and  a  is  the  fiber  direction  as  shown  in  Fig.  5.1. 
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Similar  to  the  analysis  in  the  previous  chapters,  dl5  d2  and  d3  can  be  identified  from  comparing 
the  strain  energy  densities  based  on  peridynamic  analysis  and  those  based  on  classical 
mechanics. 

Consider  a  composite  plate  with  the  fibers  oriented  in  a° 
strain  field 

^XX  — 

£yy  £"l 

Yxy  =  Y 12 

The  three  components  are  independent  of  one  another. 

For  a  bond  in  9  direction,  and  connected  to  a  point  x  in  the  domain,  the  bond  force  should  be 

/  =  c  ■  (£x  cos  9 2  +  e2  sin  9 2  +  y12  sin  9  cos  9 )  (5.7) 

From  Eqn.  4.7,  the  strain  energy  in  the  bond  becomes 

wb  =fJ-  =  cs2%/2  (5.8) 

Substituting  Eqn.  5.3  and  Eqn.  5.7  into  Eqn.  5.8  and  integrating  wb  over  the  horizon,  the  strain 
energy  density  at  the  point  x  should  be 

W  =  \  /  wb  dA  =  /027r£^  d9d$  =  537r{16(d1  +  d2)(e1  -  e2)(e i  + 

f2)  cos  2a  +  di[(fi  —  £2)2  —  7i2]  cos  4a  +  2[(3d1  +  4d2  +  8d3)(3 e\  +  2 £t£2  +  3e2  +  y?2)  + 
8(d1  +  d2)(£x  +  £2)Yi2  sin  2a  +  d1(s1  —  £2)Yi2  sin4a]}/768  (5.9) 


direction  and  subjected  to  the  following 

(5.4) 

(5.5) 

(5.6) 
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Eqn.  5.9  can  be  simplified  to  find  the  coefficient  of  each  independent  term,  as  shown  in  Table 
5.1.  The  simplification  is  achieved  based  on  Mathematica  [5], 


Table  5.1  Simplified  Eqn.  5.9  in  terms  of  independent  terms 


Coefficients 

Independent  tenns 

d1537r/96 

A  2 

cos  a *  e{ 

d1537r/96 

4  2 

cos  a  £2 

— d1537r/48 

cos  a4  £2 

— d1537r/96 

cos  a4  y\ 2 

(3dx  +  4d2)537r/96 

cos  a 2  f ! 

(— 5dx  —  4d2)d37r/96 

cos  a2 

d183n/48 

cos  a2  exe2 

d183n/96 

cos  a2  Y12 

(3  d1  +  8  d2  +  48d3)537r/768 

4 

(35d!  +  40d2  +  48d3)537r/768 

4 

(5di  +  8d2  +  16d3)537r/384 

£i£2 

(5  d1  +  8  d2  +  16d3)537r/768 

V12 

(3d!  +  4d2)83n/96 

since  cosay12£1 

(5di  +  4d2)537r/96 

sin  a  cos  ay12£2 

d183n/48 

sin  a  cos  a3  yi2£i 

— d183n/48 

sin  a  cos  a3  y12£2 

A 


On  the  other  hand,  from  the  theory  of  composite  materials  [6],  the  strain  energy  density  under 


E-XX  0 1  £yy  £ 2  ‘ind  Yxy  fl2 

W  =  \oxxex  +  i CTyy£2  +  \axyy12  (5.10) 

The  stresses  in  Eqn.  5.10  can  be  calculated  by 


®xx 

Qxx  Qxy  Qxs 

E-xx 

Oyy 

= 

Qxy  Qyy  Qys 

Eyy 

Pxy . 

Qxs  Qys  Qss 

.Yxy. 

where, 

Qxx  =  Qn™4  +  Q22n4  +  2  m2n2Q12  +  4  m2n2Q66  (5.12) 

Qyy  =  Qnn4  +  Q2  2m4  +  2  m2n2Q12  +  4  m2n2Q66  (5.13) 

Qxy  =  Qn?n2n2  +  Q22m2n2  +  (m4  +  n4)Q12  -  4 m2n2Q66  (5.14) 

Qxs  —  rrv’nQ11  —  mn3Q22  —  mn(m2  —  n2)Q12  —  2  mn(m2  —  n2)Q66  (5.15) 

Qxs  —  mri3Q  ii  —  m3nQ22  +  mn(m2  —  n2)Q12  +  2  mn(m2  —  n2)Q66  (5.16) 

Qss  —  Tn^r^Qn  +  m2n2Q22  —  2m2n2Q12  +  (m2  —  n2)2Q66  (5.17) 


m  —  cos  a  (5.18) 

n  —  sin  a  (5.19) 


Qn  — 
Q22  — 
Ql2  — 


Ei 

l-U12U2i 

E2 

l—VuVn 

E2v12 

1—V12V21 


(5.20) 

(5.21) 

(5.22) 
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(?66  —  &12 


(5.23) 


Ei_,  E1,  v12  and  G12  are  the  four  material  properties  of  orthotropic  materials.  Substituting  Eqns. 
5.11-5.23  into  Eqn.  5 . 1 0,  it  yields 

w  =  7(-  2  +  Ei4E2  +  E1  G12y12  +  2 £xExe2E2vX2  -  E2G12yl2vl2 )  cos  a4  + 

2Yi2[e2(E1(2G12  +  E2(— 1  +  v12)^  —  2.E2G12yl2 )  +  fi(E2  +  2E2G12yl2  —  E1(2G12  T 
E2v12))]  cos  a3  sin  a  +  [2 e1E1e2(E1  +  E2)  4- E\y\2  +  2eIE1E2v12  +  2 E2G12yl2vl2  + 

Ei{-2 G12yl2  +  E2 (y\2  +  2 e\v12  -  2y^2i712})]  cos  a 2  sin  a2  +  2y12[E12£2  +  2(-e1  + 
e2)E2G12v12  +  E1  (fi(2G12  +  E2(- 1  +  u12))  -  £2(2G12  +  £2Vi2))]  cos  a  sin  a3  +  [£^ef  - 

E2G12Yi2^i2  +  Ex{e\E2  +  G12yl2  +  2£1£2E2v12)]  sin  a4  +  Oh  -  £2)2G12(E1  -  E2u22)  sin  2a2} 

(5.24) 

After  simplification  processes,  the  coeffiencet  of  each  independent  tenns  can  be  found.  They  are 
listed  in  Table  5.2. 


Table  5.2  Simplified  Eqn.  5.12  in  terms  of  independent  terms 


Coefficients 

Independent 

tenns 

E\  +  4E2G12vl2  +  E1(E2  —  4G12  —  2E2v12) 
2(El  -  E2v2V2) 

4  2 

cos  a 

El  +  4  E2G12v>l2  +  E1(E2  —  4  G12  —  2E2v12) 
2(El  -  E2v'l2) 

4  2 

cos  a 

E\  +  4E2Gi2vf2  +  E1(E2  —  4  G12  —  2  E2v12) 

cos  a4  exe2 

(E1  -  E2v\2 ) 
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E\  +  4E2G12yl2  +  Ei  (E2  4  G12  2E2Vy^) 

cos  a4  7i2 

2{E1  —  £2^12) 

E\(2G12  +  E2{—  1  +  V12))  —  2E2G12Ui22 

Ei  ~  E2v\2 

cos  a 2  si 

~El  +  2E1G12  +  EiE2v  12  —  2E2Gi2vf2 

Ei  ~  E2  Vy2 

cos  a 2  sf 

E 1  +  4E2Gi2vf2  +  Ei(E2  —  4G12  —  2E2v12) 

Ei  ~  E2v\2 

cos  a 2  sxs2 

E 1  +  4E2Gi2vf2  +  Ei  (Zi2  —  4G12  —  2  E2v12) 
2(£i  -  E2vj2) 

cos  a2  7i2 

EiE2 

2  (Ei  —  E2vf2) 

ei l 

El 

2  (Ei  —  E2vf2) 

el 

E1E2V12 

El  ~  E2  Vy2 

£i£2 

E1G12  —  E2G12vl2 

2  (Ei  —  E2Ui2) 

Y 12 

Ei(2G12  +  E2(-l  +  v12)')  —  2E2G12vl2 

Ei  ~  E2v\2 

sin  a  cos  ay12s1 

El  +  2E2G12yl2  +  Ei(— 2Gi2  —  E2v12) 

Ei  ~  E2vl2 

sin  a  cos  a  y12s2 

El  +  4E2Gi2vf2  +  Ei(E2  —  4Gi2  —  2E2Vi2) 
(Ei  —  E2  Vy2) 

sin  a  cos  a3  y12s1 

El  +  4E2Gi2vf2  +  Ei(E2  —  4Gi2  —  2  E2v12) 

sin  a  cos  a3  y12s2 

(Ei  —  E2Vy2) 

Eqn.  5.9  should  be  equal  to  Eqn.  5.24.  Hence,  it  is  possible  to  set  the  coefficients  of  identical 
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terms  equal  to  each  other  and  identify  dl5  d2  and  d3. 


d^  — 


24 

S3E1n 


E i  +  E-^E2  —  + 


E  i  E1E2+4G1 


+  V  E-l  E2yj  E1E2  +  4G{ 


dy  — 


S3E1n  ' 


-3  E\  —  SE1E2  +  48  E±G12 


3E±  E1E2  +  4G 


5v/£’i£’2V  ExE2  +  4Gj 


d3 


3 

2S3E1n 


( E 1  +  SE1E2  —  20£,1fj12 


El JeiE2+4Gi2  - -  - - — 

■d - - +  4Gi2 


V12 


-E1E2+^E1E2  E1E2+4G 


2Eo  G- 1 


2 

12 


(5.25) 


(5.26) 


(5.27) 

(5.28) 


The  composite  model  has  three  parameters  and  it  can  model  shearing  deformation.  The  bond 
stiffness  is  a  continuous  function  so  the  stiffness  at  any  direction  can  be  calculated.  It  can  be 
found  that  the  bond  becomes  stronger  as  the  angle  between  the  bond  and  the  fiber  becomes 
smaller.  Similar  to  the  bar  model  in  Section  4.1,  the  Poisson’s  ratio  of  this  model  is  fixed  and  is 
related  to  the  remaining  three  independent  material  properties.  For  example,  the  Poisson’s  ratio 
of  Kevlar/Epoxy  is  0.34  while  the  Poisson’s  ratio  calculated  from  this  model  is  0.39. 


5.2  Beam  model  for  orthotropic  materials 

To  accomondate  the  four  material  properties  of  orthotropic  materials,  a  beam  model  based  on 
Section  4.2  is  proposed  here. 


From  the  composite  theory,  there  are  four  independent  material  properties,  E1,E2,v12  and  C12. 
The  composite  stiffness  varies  with  fiber  orientation.  Bond  functions  are 


(5.29) 


Cl(u' 2-u' l) 


r 


(5.30) 


They  are  identifcal  to  Eqn.  4.23  and  Eqn.  4.24.  However,  cx  and  c2  for  a  bond  in  9  direction  will 
be  dependent  on  fiber  orientations  as  follows, 

cx  =  dx  cos(0  —  a)4  +  d2  cos(0  —  a)2  +  d3  (5.31) 

c2  =  d4  (5.32) 

where  a  is  the  orientation  of  the  fiber,  as  shown  in  Fig.  5.3.  The  coefficients  dx,  d2,  d3  and  d4are 
four  independent  material  properties.  They  are  related  to  the  four  material  properties  defined  in 
the  composite  theory. 

Consider  a  composite  plate  with  fiber  in  a  direction 

fix  —  ^1 

£yy  ‘C2 

Yxy  ~  Yl2 

The  three  strains  are  independent  from  each  other. 

Similar  to  Eqn.  4.29,  the  strain  energy  in  one  bond  becomes 

c ,  (si  r  cos  62+£7r  sin  d2)2  c7(-e •,  r  cos  6  sin  6+s7r  cos  9  sin  9 )2 

wb  =  - - - -  +  - r-2 - -  (5.36) 

0  2 r  2 r3  v  J 


and  subject  to  the  following  strain  field 

(5.33) 

(5.34) 

(5.35) 


Integrate  Eqn.  5.36  to  find  the  strain  energy  density  for  a  point 
W  =  \$wb  dA  =  ^nwbrdGdr  =  ^7r5{16(d1  +  d2)82(e \  -£|)cos2a  + 
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d 182(£2  —  2s1£2  +  £2  —  Yi2)  cos  4 a  +  2[24d4(s1  —  e2)2  +  9d182ef  +  12d282s2  + 

24d352£f  +  6  d182£1s2  +  8d282£1£2  +  16d352f:1£2  +  9d182£2  +  12d252£f  +  24  d382£2  + 
24d4y22  +  3d152y22  +  4d252y|2  +  8d352y22  +  8(d4  +  d2)52(£1  +  £2)Yi2  sin  2a  + 
d182(£1  —  £2)Yi2  sin  4a]}  (5.37) 


The  Strain  energy  density  based  on  the  composite  theory  is  the  same  as  Eqn.  5.24.  Similar  to 
Section  5.1,  by  setting  Eqn.  5.24  equal  to  Eqn.  5.37  and  comparing  the  coefficients  of  the 
independent  terms,  the  following  equations  are  obtained 

48  (E4  +E1E2—4:E1G12  —  2E1E2V12  +  4E2G12V^2) 


dy  = 


n(E1-E2vl2)83 


.  12  (3 E3  4  SE4E2  —  16&L  G4 2 — 8£4£,2V42  4 1 GE2  ^12^12) 

2  n(E1-E2vl2)83 


do  — 


3(£’1  +  5E4E2  —  8E1G12  —  2£4£’2rli24'8£'2(ji21;’i2) 


dA  =  - 


n  {E1-E2vl2)S3 


4(— Ex  G12  ^E1E2v12 +£’2^12^12) 


n(iE1-E2vl2)S 


(5.38) 

(5.39) 

(5.40) 

(5.41) 


5.3  Calculation  of  stresses  from  peridynamics 

Similar  to  Section  4.2,  stresses  can  be  defined  and  calculated  from  peridyamics.  They  can  be 
used  for  some  special  comparison  but  not  necessary  in  peridynamic  simulations. 


Consider  a  case  with  £xx  —  £4  and  £yy  —  £2.  The  stress  oxx  can  be  calcualted  by  using  Eqn.  4.39 


and  is  given  below 


@xx 


rS  rS  rCOS  ^ 

_  J 0  't  J-cos-1- 


(fx'  ■  cos  6  —  fy'  ■  sin  0)r  dddrdt  —  —  n8{48d4(£1  —  £2)  + 
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2(3dx  +  4 d2  +  8d3)(3£1  +  £2)52  +  82[16£1  cos  2 a  (dx  +  d2)  +  d1(£1  —  £2)  cos  4a]}  (5.42) 


Substituting  Eqns.  5.38  -  5.41  into  Eqn.  5.42,  it  yields 

®xx 

2  7 {3fi E2  +  £2E2  +  3 s1E1E2  +  £2E1E2  +  4£1E1G12  —  4£2E1G12  +  2e1E1E2v12  + 

3(h1-h2V12)  '■ 

6£2E1E2u12  -  4s1E2G12u22  +  A£2E2G12vl2  +  4£1£1(£1  -  E2)  cos  2a  +  (fy  -  £2)(£2  + 
4E2G12u22  +  E1(E2  —  4 G12  —  2£'2u12))  cos  4a]  (5.43) 

While  from  the  composite  theory,  a**  can  be  expressed  as 

1  o 

0**  =  QxxZ  1  +  — ^-{EiCgiEi  +  £2^2^12)  cos  a4  +  +  4(-£i  + 

(.£l~£2v12  ) 

£2)E2G12vl2  +  E^E-2  +  4£1G12  -  4 £2G12  +  2£1E2U12)]  cos  a2  sin  a2  +  + 

£2v12)sin  a4}  (5.44) 

With  further  simplification,  it  can  be  found  that  Eqn.  5.43  is  identical  to  Eqn.  5.44. 

5.4  Laminated  plate  under  static  loading 

In  this  section,  it  is  to  verify  the  proposed  peridynamic  model  with  an  anlytical  solution.  A 
simple  tensile  test  is  perfonned  on  a  laminated  plate.  The  peridynamic  results  will  be  compared 
with  the  results  obtained  from  the  composite  theory. 

Consider  a  100  mm  x  100  mm  laminated  plate  with  fibers  in  a  direction  as  shown  in  Fig.  5.4.  A 
tensile  pressure  of  10  MPa  is  applied  at  the  bottom  and  the  top  of  the  plate.  The  plate  is  made  of 
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E-Glass/Epoxy  and  the  material  properties  are  shonw  in  Table  5.3. 


Table  5.3  Material  propteties  of  E-Glass/Epoxy 


Longitudinal  Young’s  modulus,  E1 

41  GPa 

Transverse  Young’s  modulus,  E2 

10.4  GPa 

Poisson’s  ratio,  v12 

0.28 

Shear  modulus,  G12 

4.3  GPa 

Mass  density,  p 

1970  kg/m3 

Based  on  the  composite  theory  [6],  the  components  of  the  complicance  matrix  are 


Sn  -  1  /£i 

(5.45) 

“-*22  =  1/^2 

(5.46) 

‘-*12  ~  ~V12/El 

(5.47) 

$66  ~  1/^12 

(5.48) 

Strains  from  the  composite  theory  can  be  calcuated  as  follows 


^xx 

r  s  s  si 

°xy  °xs 

®xx 

£yy 

= 

s  s  s 

Jxy  ^yy  uys 

Gyy 

_Yxy. 

s  s  s 

yxs  °ys  lJss_ 

Pxy . 

(5.49) 


The  transformed  compliance  matrix  is  caculated  as 


s 

°xx 

s 

°xy 

s  1 

^xs 

Sn 

‘-’12 

0  - 

s 

°xy 

s 

s 

0 ys 

=  t~1 

^12 

$22 

0 

-S 

[2  xs 

-S 

2  ys 

—  s 

2  ss\ 

0 

0 

-S 

2^66] 

(5.50) 


where 
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T  = 


2 

771 

n2 

2  mn 

n2 

m2 

—2  mn 

-—mn 

mn 

m2  —  n‘ 

and  m  =  cos  a  and  n  =  sin  a. 


(5.51) 


From  Eqn.  5.49,  the  displacement  field  can  be  obtained  from  the  composite  theory. 
Displacements  from  peri  dynamics  are  compared  with  those  from  the  composite  theory  in  Fig. 
5.5-Fig.  5.10  with  a  —  0°,  45°  and  60°.  As  can  be  seen,  the  results  from  peridynamics  and  those 
from  composite  theory  are  identical  to  each  other. 

5.5  Free  vibration  of  a  laminated  beam 

The  free  vibration  of  a  laminated  beam  is  investigated  in  this  section  by  the  Classical  Laminated 
Beam  Theory  and  the  solution  will  be  used  to  verify  that  obtained  from  peridynamic  model. 

Consider  a  simply  supported  beam  as  shown  Fig.  5.11.  The  length  of  the  beam  is  L  =  a  ■  h  and 
the  thickness  of  the  beam  is  h  —  10  mm,  where  a  is  the  aspect  ratio  of  the  beam.  The  beam  is 
made  of  Kevlar/Epoxy  with  fibers  oriented  in  x  direction.  The  material  properties  are  shonw  in 
Table  5.4. 


Table  5.4  Material  propteties  of  Kevlar/Epoxy 


Longitudinal  Young’s  modulus,  E1 

80  GPa 

Transverse  Young’s  modulus,  E2 

5.5  GPa 

Poisson’s  ratio,  v12 

0.34 

Shear  modulus,  G12 

2.2  GPa 

13 


Mass  density,  p 


1380  kg/m 3 


From  [7],  the  governing  equations  of  the  beam  are 

.  d2u°  „  d3w  ?  n  „ 

B ^  P“2u°  =  ° 

(5.52) 

d4w  d3u°  2  n 

P°1W  =  0 

(5.53) 

where 

—  /-ft/2  =  Qnh 

(5.54) 

Bn  =  (?n  zdz  =  0 

(5.55) 

£ii  =  I%22CiiZ2  dz  =  Q nh3/12 

(5.56) 

A  solution  satisfying  the  governing  equations  is 

u0  —  U  cos  px 

(5.57) 

w  =  I/F  sinpx 

(5.58) 

Eqn.  5.57  and  Eqn.  5.58  satisfy  the  simply  supported  boundary  conditions  automatically. 
Substituting  Eqn.  5.57  and  Eqn.  5.58  into  Eqn.  5.52  and  Eqn.  5.53,  it  yields 

-Atlp2U  -  pco2U  =  0  (5.59) 

Dnp4lE  -  pai2W  =  0  (5.60) 


Solving  Eqn.  5.59  and  Eqn.  5.60,  it  can  be  concluded  that 


co 


2 


QnP4h3 
12  p 


(5.61) 
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U  =  0 


(5.62) 


If  the  initial  condition  of  the  beam  is 

w(t  =  0)  =  1  x  10_5sin^x  (5.63) 

then 

W  =  lx  1(T5  (5.64) 

P  =  \  (5-65) 

The  solution  of  the  problem  should  have  the  following  forms 

w'(x,  t)  =  W  sin  px  cos  cot  (5.66) 

u'(x,  z,  t)  =  —  zpVK  cos  px  cos  cot  (5.67) 

where  W,  p  and  co  can  be  found  from  Eqn.  5.64,  Eqn.  5.65  and  Eqn.  5.61,  respectively. 

The  same  problem  can  be  simulated  by  peridynamics.  The  displacement  history  of  point  A  and 
point  B  (Fig.  5.11)  are  recorded.  Fig.  5.12  compares  the  peridynamic  results  with  those  from  the 
composite  beam  theory  for  the  aspect  ratio  a  =  5.  Fig.  5.13  compares  the  peridynamic  results 
with  those  from  the  composite  beam  theory  for  the  aspect  ratio  a  =  20.  Results  from  the  two 
methods  show  good  match  in  vertical  displacement  w  of  point  B.  For  horizontal  displacement  it 
of  point  A,  peridynamic  result  is  almost  the  same  as  the  beam  theory  result  when  the  aspect  ratio 
a  —  20.  However,  there  is  difference  between  the  peridynamic  result  and  the  beam  theory  result 
when  the  aspect  ratio  a  —  5.  This  is  because  the  beam  theory  assumes  no  variation  of  vertical 
displacement  when  the  beam  is  slender.  With  a  small  aspect  ratio,  such  as  a  —  5,  this  variation  is 
not  neglible  and  the  beam  theory  does  not  provide  a  good  approximation. 
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5.6  Comparisons  of  crack  propagation  velocity  with  experiments 

Experimental  studies  on  dynamic  crack  propagation  in  fiber-reinforced  composite  materials  have 
been  conducted  by  Zheng[9],  Rosakis[8],  Stout[10]  and  Coker[ll,12].  It  has  been  shown  that  a 
weak  fracture  plane  usually  occurs  between  fiber  and  matrix  in  unidirectional  fiber-reinforced 
composites.  Due  to  material  anisotropy,  the  wave  speed  along  the  fiber  direction  is  very  different 
from  that  along  the  perpendicular  direction.  Dynamic  crack  propagation  has  been  commonly 
investigated  by  finite  element  method.  A  limited  number  of  computational  studies  have  been 
reported  by  Huang  [13],  Hwang  [14],  Kumar  [15],  Stout  [10],  Lo  [16],  Sun  [17]  and  Pandey[18]. 
The  limit  of  computational  works  is  likely  due  to  the  requirement  of  remeshing  and  the 
complexity  of  handling  elements  once  crack  starts.  In  this  section,  dynamic  crack  propagation  in 
a  unidirectional  graphite/epoxy  composite  is  studied  with  the  use  of  peri  dynamics.  The 
computational  results  are  validated  by  the  experimental  results  given  in  refefrence  [8]. 

Consider  a  76  mm  x  152  mm  unidirectional  graphite/epoxy  fiber-reinforced  composite  plates 
under  three-point  bending  [8]  as  shown  in  Fig.  5.14.  The  fiber  is  in  0°  direction.  The  material 
properties  are  shown  in  Table  5.5  [8,  22].  There  is  a  notch  with  a  length  of  15.2  mm,  i.e.  20%  of 
the  plate  width,  at  the  left  boundary  of  the  plate.  This  crack  length  is  used  because  it  was  used  in 
the  past  to  produce  reliable  results  in  dynamic  fracture  experiements  [19].  To  minimize  residual 
stresses  due  to  machining,  a  low-speed  diamond  saw  was  used  to  produce  the  initial  notch  with  a 
widtch  of  approximately  1.5  mm. 

Table  5.5  Material  propteties  of  graphite/epoxy 
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Fongitudinal  Young’s  modulus,  E1 

150  GPa 

Transverse  Young’s  modulus,  E2 

11.6  GPa 

Poisson’s  ratio,  v12 

0.36 

Shear  modulus,  G12 

3.5  GPa 

Mode  I  intralaminar  fracture  energy  for  longitudinal  loading,  G10 

77.9  J/m2 

Mode  I  intralaminar  fracture  energy  for  transver  loading,  G2Q 

5  J/m2 

Mass  density,  p 

1590  kg/m3 

A  drop-weight  tower  is  used  to  introduce  impact  on  the  opposite  side  of  the  notch  with  an 
impacting  speed  of  v0  —  4  m/ s.  After  the  impact,  stress  waves  propagate  to  the  interior  of  the 
plate  and  then  reflects  from  the  boundaries.  Because  of  the  anisotropy  of  the  material,  stress 
waves  travel  in  different  directions  at  different  velocities.  Experiments  show  that  the  crack  starts 
to  propagate  at  about  25 ps  after  the  impact  so  the  effects  of  dispersion  are  not  very  important 
since  the  applied  stress  pulse  is  very  long  (about  120  ps)  compared  with  the  time  of  crack 
initiation.  Therefore,  loading  is  continuously  applied  throughout  the  entire  event. 

The  real-time  visualization  of  dynamic  fracture  is  produced  by  an  optical  method  of  coherent 
gradient  sensing  (CGS)  in  reflection  [20,  21].  Imaging  is  perfonned  with  a  rotating-mirror  high¬ 
speed  camera.  Details  of  the  CGS  system  can  be  found  in  [8,  20,  21]. 

Fig.  5.15  shows  the  crack  propgation  velocity  from  the  peridynamic  simulation.  The  initial  ime 
t  =  0  is  used  to  denote  the  beginning  of  the  crack  propagation.  For  negtive  time,  v  =  0  m/s. 
The  crack  starts  from  about  700  m/s  and  accelerates  to  900  m/s  within  about  10  fis.  It  then 
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decelerates  to  less  than  500  m/s  in  40  /rs  after  the  initiation.  This  deceleration  is  believed  to  be 
due  to  the  fact  that  the  growing  crack  tip  enters  a  region  of  high  compressive  stress  zone  as  it 
appoaches  the  loading  area.  The  peridynamic  computational  results  are  compared  with  the 
experimental  results  from  Ref.  [8].  As  shown  in  Fig.  5.15,  they  agree  each  other  reasonably  well. 

5.7  Dynamic  fracture  mode  in  unidirectional  composites 

In  order  to  investigate  the  behavior  of  cracks,  Wu  [24]  conducted  experiments  with 
unidirectional,  fiberglass-reinforced  Scotch  composites  with  a  centered  precrack  in  the  direction 
of  fibers.  The  composites  were  loaded  with  tension,  pure  shear  and  combined  tension  and  shear. 
In  all  three  cases,  it  was  observed  that  the  crack  propagated  in  the  same  direction  as  the  fiber 
direction.  Finite  element  analysis  were  also  used  to  study  the  damage  path  and  failure  initiation 
of  pre -notched  composites  by  Boger  [25]  and  Satyanarayana  [26].  They  prediected  damge  in 
composite  plates  notched  in  the  center  for  different  layups  under  tension.  Both  experimental 
results  and  simulation  results  showed  that  the  crack  path  and  failure  initiation  depends  on  fiber 
orientation. 

In  this  section,  the  crack  propagation  path  and  dynamic  fracture  mode  of  unidirectional  fiber- 
reinforced  composites  are  studied  using  the  proposed  peridynamic  model.  The  qualitive 
comparison  of  the  peridynamic  results  with  those  from  experiments  are  of  interest. 

Consider  the  compact  tension  test  on  a  100  mm  x  200  mm  carbon/epoxy  unidirectional 
composite  plate  with  a  20  mm  pre-notch  at  the  center  as  shown  in  Fig.  5.16.  The  plate  is  loaded 
at  the  top  and  the  bottom  boundary  by  a  uniform  stress  er.  The  fiber  orientation  is  a.  The  material 
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properties  of  the  carbon/expoxy  plate  are  shown  in  Table  5.6  [23]. 


Table  5.6  Material  propteties  of  graphite/epoxy 


Longitudinal  Young’s  modulus,  E1 

329  GPa 

Transverse  Young’s  modulus,  E2 

6  GPa 

Poisson’s  ratio,  v12 

0.346 

Shear  modulus,  G12 

4.4  GPa 

Mode  I  intralamina  fracture  energy  for  longitudinal  loading,  G10 

15.49  k  J/m2 

Mode  I  intralamina  fracture  energy  for  transver  loading,  G2o 

0.168  kj/m2 

Mass  density,  p 

1630  kg/m3 

Fig.  5.17,  Fig.  5.18,  Fig.  5.19,  Fig.  5.20  and  Fig.  5.21  show  the  peridynamic  simulation  results 
for  a  =  0°,  a  =  30°,  a  =  45°,  a  =  60°  and  a  =  90°,  respectively.  In  each  figure,  there  are  three 
contour  plots  with  (a)  vertical  diaplacement  of  the  plate,  (b)  strain  energy  density  distribution  and 
(c)  local  damge  of  the  plate.  Different  color  bars  are  associated  with  different  plots  with  red 
indicating  the  highest  value  while  blue  the  lowest  value.  Crack  paths  can  be  clearly  seen  from  the 
strain  energy  density  contour  plot  as  there  are  always  strain  energy  concentraions  at  the  crack 
tips.  The  local  damage  is  defined  by  Eqn.  2.12  and  Enq.  2.13,  which  show  different  levels  of 
damage.  A  proper  cuttoff  value  can  be  defined  to  judge  if  there  is  a  crack.  In  all  cases,  the  crack 
propagates  in  the  same  direction  as  the  fibers,  which  is  consistent  with  the  experimental 
observations  from  Wu  [24].  The  damage  is  due  to  the  seperation  between  matrix  and  fiber.  There 
is  no  fiber  breakage. 
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As  expected,  in  the  a  —  0°  case,  the  crack  propagates  in  the  same  direction  as  the  pre-notch. 
This  also  matches  with  the  computational  and  experimental  results  in  Section  5.6.  In  the  smaller 
angle  case  a  —  30°,  aside  from  the  major  crack,  which  propagates  in  30°,  there  is  matrix 
shattering  at  the  sides  of  the  plate  in  0°  direction.  The  matrix  shattering  happens  before  the  crack 
starts  to  propagate  as  shown  in  Fig.  5.22.  It  starts  at  the  lateral  of  the  plate  and  propagates  to  the 
interior  of  the  plate.  From  Fig.  5.18  (c),  the  matrix  shattering  is  not  as  severe  as  the  major  crack 
since  only  20%  of  the  bonds  are  broken.  However,  in  the  major  carck,  more  than  70%  of  the 
bonds  are  broken.  This  is  why  matrix  shattering  is  only  a  material  softening  and  may  not  be  seen 
from  experimental  observation  as  reported  in  [24],  For  the  a  —  90°  case,  the  composite  plate 
fails  due  to  splitting  caused  by  shear  stress  in  the  matrix.  It  which  matches  with  the  findings  of 
Boger  [25]. 

5.8  Conclusions 

A  peridynamic  orthotropic  material  model  baesd  on  the  beam  model  is  proposed  in  this  chapter. 
There  are  four  independent  mateiral  parameters  in  this  model  and  it  matches  with  the  four 
material  properties  for  two-dimensional  orthotropic  materials.  The  bond  material  properties 
depend  on  these  four  mateiral  parameters  and  the  angle  between  the  bond  orientation  and  fiber 
orientation.  This  results  in  the  continuity  of  the  bond  stiffness  function  with  no  need  of 
remeshing  for  different  fiber  orientations.  This  model  is  verified  by  a  static  tensile  test  and  a 
vibration  problem  of  a  laminated  beam. 

Dynamic  damage  propagation  problems  in  composite  materials  can  br  greatly  benefited  from 
peridynamics.  The  prediction  of  damage  initiation  and  crack  proagation  of  composite  materials  is 
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complex  using  traditional  methods,  such  as  finite  element  analysis,  due  to  its  anisotropy.  As 
investigated  in  peridynamic  simulations,  there  is  no  need  of  tracking  each  crack  propagation, 
finding  different  damage  modes  and  applying  different  damage  rules.  Damage  happes 
automatically.  A  single  edge  notch  test  is  simulated  in  this  chapter  and  the  results  match  with  the 
experimental  results.  Crack  path  and  failure  initiation  of  a  center  notch  plate  is  predicted 
successfully  by  peridynamics  when  comparing  with  the  experimental  results. 
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Figure  5. 1  A  composite  plate  with  fiber  in  a°  direction  and  a  bond  in  9°  direction. 
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Figure  5.3  Beam  model  for  orthotropic  materials. 
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Figure  5.4  Pulling  test  in  a  composite  laminate. 
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Figure  5.5  Comparison  of  ux  of  0°  laminate  calculated  from  peridynamics  (top)  and  composite 
theory  (bottom). 
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Figure  5.6  Comparison  of  ux  of  45°  laminate  calculated  from  peridynamics  (top)  and  the 
composite  theory  (bottom). 
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Figure  5.7  Comparison  of  ux  of  60°  laminate  calculated  from  peridynamics  (top)  and  the 
composite  theory  (bottom). 
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Figure  5.8  Comparison  of  uy  of  0°  laminate  calculated  from  peridynamics  (top)  and  the 
composite  theory  (bottom). 
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Figure  5.9  Comparison  of  uy  of  45°  laminate  calculated  from  peridynamics  (top)  and  the 
composite  theory  (bottom). 


32 


0  0.02  0.04  0.06  0.08  0.1 


0  0.02  0.04  0.06  0.08  0.1 


0.08- 


0.06- 


0. 04- 


0.02- 


Figure  5.10  Comparison  of  uy  of  60°  laminate  calculated  from  peridynamics  (top)  and  the 
composite  theory  (bottom). 
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Figure  5.11  Free  vibration  of  a  laminated  beam  with  fibers  in  x  direction. 
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Figure  5.12  Peridynamic  results  compared  with  composite  beam  theory  results  with  a  —  5. 
Top:  horizontal  displacement  it  of  point  A.  Bottom:  vertical  displacement  w  of  point  B. 
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Figure  5.13  Peridynamic  results  compared  with  composite  beam  theory  results  with  a  —  20. 
Top:  horizontal  displacement  it  of  point  A.  Bottom:  vertical  displacement  w  of  point  B. 
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Figure  5.14  An  unidirectional  composite  plate  with  single  edge  notch  under  three  point  bending. 
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Figure  5.15  Comparison  of  crack  propagation  velocity  between  peri  dynamics  and  experiment. 
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Figure  5.16  Compact  tension  test  for  a  unidirectional  composite  plate. 
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Figure  5.17  Simulation  results  for  t  —  50  fis 
and  a  —  0°: 

(a)  vertical  displacement,  (b)  strain  energy 
density,  (c)  local  damage. 
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Figure  5.18  Simulation  results  for  t  =  70  /rs 
and  a  —  30°: 

(a)  vertical  displacement,  (b)  strain  energy 
density,  (c)  local  damage. 
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Figure  5.19  Simulation  results  for  t  =  70  [is 


and  a  —  45°: 


(a)  vertical  displacement,  (b)  strain  energy 
density,  (c)  local  damage. 


Figure  5.20  Simulation  results  for  t  =  90  [is 
and  a  —  60°: 

(a)  vertical  displacement,  (b)  strain  energy 
density,  (c)  local  damage. 
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Figure  5.21  Simulation  results  for  t  = 
100.5  /rs  and  a  —  90°: 

(a)  vertical  displacement,  (b)  strain  energy 
density,  (c)  local  damage. 
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Figure  5.22  Simulation  results  for  t  —  40  /rs 
and  a  —  30°: 

(a)  vertical  displacement,  (b)  strain  energy 
density,  (c)  local  damage. 
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